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ABSTRACT 

Supermassive black holes in the centres of giant elliptical galaxies are thought to 
be capable of inducing chaos and eliminating or preventing triaxiality in their hosts if 
they are sufficiently massive. We address whether this process operates in real systems, 
by modeling the stellar kinematics of the old elliptical NGC 4365. This galaxy has a 
mean stellar population age > 12 Gyr and is known for its kinematically decoupled 
core and skew rotation at larger radii. We fit the two-dimensional mean velocity field 
obtained by the SAURON integral-field spectrograph, and the isophotal ellipticity 
and position-angle profiles, using the velocity field (VF) fitting approach. The models 
constrain the system's intrinsic shape between 0.03 and 0.5 effective radii, as well as 
its orientation in space. We find NGC 4365 to be strongly triaxial ((T) « 0.45) and 
somewhat flatter than it appears {{c/a) « 0.6). Axisymmetry or near axisymmetry 
(T < 0.1) is ruled out at > 95% confidence for 1'.'6 < R< 3'.'2 (0.03 < R/Ve < 0.06), 
and at > 99% confidence at larger radii. There is an indication of an outward triax- 
iality gradient. The line of sight is constrained to two narrow bands on the viewing 
hemisphere. In the most probable orientation the long axis points roughly toward the 
observer, extending to the southwest in projection. The stellar population age implies 
that strong triaxiality has persisted for hundreds of dynamical times. This rules out 
black holes > 3 x 10^ Mq, which numerical simulations indicate would either have 
globally axisymmetrized the galaxy or made the inner several arcseconds spherical. 
The Mbh-ct relation predicts A/bh ~ 4 x 10^ Mq, which would probably not preclude 
long-lived triaxiality and is consistent with the observations. There must also be an 
unequal population of direct and retrograde long-axis tube orbits outside the kinemat- 
ically decoupled core. This, combined with the small isophotal twist, limits the rate of 
figure rotation (tumbling) about the short axis, and places corotation at > 8 effective 
radii. NGC 4365 lends support to a picture in which supermassive black holes, though 
omnipresent in luminous giant elliptical galaxies, are not massive enough to alter their 
global structure through chaos. 

Key words: black hole physics — stellar dynamics — galaxies: elliptical and lenticular, 
cD — galaxies: evolution — galaxies: individual (NGC 4365) — galaxies: kinematics and 
dynamics 



1 INTRODUCTION 

There is growing recognition that the dynamical structure 
of hot stellar systems — elliptical galaxies and the bulges 
of spirals — is intimately connected with the supermassive 
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black holes that reside at their centres. The most com- 
pelling evidence for this connection is the tight correlation 
between black hole masses and the velocity dispersions of 
their host galaxi es measured far from the black hole. The 
Mbh -t relation llFerrarese fc Merrittl|2000l : iGebhardt et alJ 
l200Ci) indicates that central black holes somehow "know" 
about the properties of their hosts, and has provoked specu- 
lation t hat host dynamics control the growth of the central 
object JEI Zant et al.ll2003l : lA'dams et al.ll2003ll . But the ex- 
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act mechanisms by which this control might be exerted are 
still obscure. 

By contrast, a mechanism by which a central black 
hole may, in a turnabout move, affect the structure of its 
host IS understood, or at least identified. Important fam- 
ilies of regular orbits that pass close to the centre of the 
system ca n be rendered chaotic by scattering off the cen- 
tral mass ('L ake fc NormanlllQsA iGerhard fc Binnevlll985l: 
iGerhard 198a) . It has thus been suggested that growth of 
a black hole can lead to the destruction of bars in spirals 
JHasan fc Nor man 1 99(|; iNorman et aljl 996|) or to the elim- 
ination of tria xiality in ellipticals iValluri fc Merrittill998l : 
iMerritt fc QuinlarJl998l) . In this case, however, what is lack- 
ing is direct evidence that this mechanism has played a role 
in real systems. 

When central point masses or density cusps are added 
to triaxial potentials, box orbits are transformed to st ochas- 
tic orbits or resonant "boxlets" llSchwarzschildlll993r .^ De- 
spite some speculation that the l oss of boxes would pre- 
clude nearly all triaxial equilibria llMerritt|ll99dl . it is now 
understood that stochastic orbits can contribute construc- 
tively to maintaining the triaxial figure (jSchwarzschild 1 9931 : 
iMerritt fc Fridmanlll99a) . The importan t questions now fo- 
cus o n the chaotic mixing time scale iSiopis fc KandrutJ 
I2OOC1) . and whether, for practical purposes, one can con- 
sider there to be many distinct, slowly diffusing, stochas- 
tic orbits at each energy, rather than only one such or- 
bit o ccupying the entire Arnol'd web (IMerritt fc FridmanI 
| l99d^. Studies of indi vidual orbits (IVaUuri fc Merrittill998l: 
IPoon fc Merritt P2001ll as well as evolutionary A'^- body simu- 
lations ^Merr itt fc Quin lan""l998': ISellwood> "2002') agree that 
the box-orbit region of phase space becomes fully chaotic 
when the central point mass reaches ~ 1% of the total 
system mass. This transition happens somewhat earlier for 
more elongated figures or steeper cusps. At lower masses, a 
fully chaotic zone appears over a limited range of energies, 
affecting a region whose enclosed s tellar mass is of order 
10 times that of the central o bject llPoon fc Merritti 120021 : 
iHollev-Bockelmann et 811120021) . 

The current Mbh-c relation (|TYemain^^^aL|_|2003i, 
coup l ed with the Fundam ental Plane dPiorgovski fc Davis! 
ll987l : lDressler et aLlll987^ ■ implies that black holes wiU typ- 
ically fall about a factor of 10 short of the 1% mass needed 
to induce global chaos and a major restructuring of the host 
galaxy. If this is universally the case, triaxiality should per- 
sist over long times, and we should expect to find ma ny old, 
triaxial ellipticals. However, IValluri fc Merritg (|l993) argue 
that the threshold for chaos should be proportionally lower 
in lower luminosity galaxies, which have steeper cusps and 
shorter dynamical times, and suggest that galaxies fainter 
than around Mb ~ —20 should have been axisymmetrized 
by chaos. Consequently, determining the actual abundance 
of triaxial systems as a function of luminosity will be an im- 
portant check on our understanding of the masses of black 
holes and their influence on their hosts. 

Determining the true shapes of elliptical galaxies cannot 
be done by photometr y alone. K inematic data are essential 
llBinnevlll985. : .Franx ct ahlliggif) , and the "traditional" ma- 



^ For some special counter-examples see lSridhaj fc Toumal Jl997l . 
l200(t. 



jor and minor axis kinematic profiles are not enough. IStatlerl 
(I1994a|) advocated a minimum of 4 long-slit position angles 
for a credible constraint on triaxiality, a requirement that is 
prohibitively costly in telescope time. Happily, this has be- 
come a moot point with the advent of integral-field spe ctro- 
graphs. The SAURON instrument dBacon et al.ll200J) was 
developed specifically to map the kinematic and spectral 
properties of galaxies, with full two-dimensional coverage 
on the sky. SAURON has been used to observe a sample of 
72 E, SO, and Sa galaxies (Ide Zeeuw et alJl2002h . producing 
an extensive catalog of kinematic and line-index measure- 
ments. In this paper we fit triaxial models to the SAURON 
kinematic dataset for the old elliptical NGC 4365, in order 
to estimate its intrinsic shape and constrain its orientation 
relative to the line of sight. 

The properties of NGC 4365 have recently been sum- 
marized by Davics ct al. (200l|). It is a luminous giant ellip- 
tical {Mb = —21.05 for a distance modulus of 31.55) with a 
very shallow central cusp. The negative logarith mic slope of 
the i nner surface brightness profile is 7 « 0.13 IIRest et alJ 
bOOlT) . Photometrically, the galaxy has near-constant ellip- 
ticity e ~ 0.26, and a modest isophotal twist of about 10° 
from 0.02 to 2 effective radii (re = 57"). The isophotes are 
discy in the inner 4" and boxy farther out. NGC 4365 is 
probably best known for its kinematics, having both a kine- 
matically decoupled core and skew rotation at larger radii. 
The apparent rotation is about the photometric minor axis 
in the core, and ~ 40° from the major axis at 20" from 
the centre. These kinematic signatures were recognized by 
ISurma fc Benden (|l99a) as hal lmarks of triaxiality . From 
the SAURON hue index maps, iDavies et all (120011) deter- 
mine that the stellar population in NGC 4365 has a mean 
age of 14Gyr, both in and outside the kinematically dis- 
tinct core, with at most a few percent admixture of an 
intermediate-age population in the innermost l"6. They find 
that putting 6% of the mass in a 5-Gyr-old population with 
the same metallicity as the rest of the galaxy ca n account for 
the sl ightly elevated central H/3 line-strength. iDavies et alJ 
d200ir) conclude that the galaxy has been largely unchanged 
in the last 12 Gyr. However, there is also evidence of a 
popu lation of intermediate-age (2-8 Gyr o ld) globular clus- 
ters llPuzia et al.ll200a iLarsen et al.ll2003^ . suggesting that 
a minor merger event a few billion years in the past may 
have modified the structure of the galaxy and possibly con- 
tributed to the kinematically decoupled core. 

In § 121 below we describe the methods we use to con- 
strain the shape and orientation of the galaxy, including 
the nature of the models, the handling of the observational 
data, and the statistical formalism. § |3 presents the basic 
results. We find that NGC 4365 is highly triaxial, with ax- 
isymmetry ruled out at better than 99% confidence between 
approximately 2" and 30". We find the line of sight to be 
constrained to two fairly narrow strips on the hemisphere 
of viewing angles. We are also able to constrain, to a lim- 
ited degree, the contributions of the main tube orbit families 
to the mean velocity field. In § 21 we show that the persis- 
tent triaxiality of this old system rules out black holes larger 
than a few 10^ Mq , but is probably consistent with predic- 
tions of the Mbh-o" relation. We also use the predominance 
of long-axis tubes in the outer part of the system and the 
lack of a strong photometric twist to limit the rate of figure 
rotation. We briefiy discuss the implications for orbit-based 
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Figure 1. Construction of a velocity field on one fiducial shell, (a) A confocal coordinate system linked to the model triaxiality (in 
this example T = 0.5) represents the streamlines of short-axis tubes (red dotted lines) and long-axis tubes (blue dotted lines). Given a 
density distribution, specifying the velocity crossing an arc in the [x, z) plane {heavy colored line) then determines the mean flow field 
everywhere [black vectors and streamlines). The boundary arc is divided into two segments (red, blue) describing the contribution of 
short-axis and long-axis tubes, respectively, (b) The function f*(t), which determines the velocity along the boundary arc, showing the 
definition of the constants C, kx, and kz. The example in (a) has (7 = 1, kx = kz = 0. 



(Schwarzschild) and moment-based (Jeans) models, before 
summarizing the principal results in § |^ 



2 METHOD 

We determine the intrinsic shape of the galaxy by fitting 
models to the mean line-of-sight velocities and the radial 
profiles of isophotal ellipticity and major-axis position an- 
gle. We do not make use of the dispersion or the higher mo- 
ments (except for a v/a constraint described below) , because 
our goal is not to determine the full three-dimensional mass 
profile. That is a job better left to the more sophisticated — 
but far slower — Schwarzschild method. Our goal is only to 
constrain the shapes of the isodensity surfaces and the ori- 
entation of the galaxy relative to the line of sight. We can 
extract this information from the mean rotation and pho- 
tometry because the mean velocities in a non-tumbling tri- 
axial galaxy arise from the populations of stars on short-axis 
(Z) and long-axis (X) tube orbits, and the geometry of these 
orbits is determined almost entirely by the triaxiality of the 
potential. 



2.1 Models 

The velocity field ( VF) fitting approach has been d escribed 
in previous papers llStatlej200ll : IStatler et alll999l . and ref- 
erences therein) . The models are based on analytic solutions 
to the equation of continuity. Though the models are approx- 
imate, their calculation is extremely fast, making it possible 
to explore a far larger parameter space than can be covered 
by more computationally intensive orbit-based methods. We 
sketch the essentials here, and explain certain details in the 
Appendix. 

The figure of the galaxy is assumed to be stationary, and 
the principal axes of the isodensity surfaces are assumed to 
be aligned throughout. For each projected radius at which 



there are data to be fitted, we calculate an internal (3- 
dimensional) mean VF by solving the equation of continuity 
on a deprojected shell. This solution involves the assumed 
density figure of the model as well as the geometry of the 
streamlines of the mean flow. To approximate the latter, we 
assume that the mean motions (not the single-particle or- 
bital velocities) of X and Z tubes follow coordinate lines in 
a confocal coordinate system that is keyed to the triaxial- 
ity T of the mass distribution^ at the radius of the shell. 
We have verifled the accuracy of this approximation for a 
wide variety of tube orbits numerically integrated in loga- 
rithmic potentials of constant T (! And erso n & Statlcr 199q). 
For potentials with large triaxiality gradients, linking the 
streamlines to the local T is probably inaccurate, but this is 
moot here since we find no evidence for a large gradient in 
NGC 4365. 

Figure ^ shows the geometry of the streamlines and 
the internal velocity field, drawn on a spherical shell in one 
model. Triaxial systems can also have nonzero radial mean 
motions, due to the elongation of tube orbits contrary to the 
elongation of the potential. We compute each velocity field 
both on a spherical shell and on an ellipsoidal shell intended 
to overestimate the ex pected r adial m otions, and average 
the results [see § 2.4 of lStatleJ l)l994Ef) for details]. 

The solution of the continuity equation also requires a 
boundary condition, which is specified at each radius as the 
mean velocity along an arc in the {x, z) plane (Fig.^). This 
boundary splits into two segments, the upper determining 
the contribution from X tubes, and the lower the contribu- 
tion from Z tubes. Specifying the velocity on the boundary 
can be seen as equivalent to specifying the orbit populations 
in the model — or more precisely, that part of the orbit popu- 
lation that determines the observable velocity field. We write 
the absolute magnitude of the boundary velocity as v*{t), 



^ The triaxiality parameter T = (a^ — b^)/(a^ 
and c are the long, middle, and short axes. 
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Figure 2. (a) SAURON mean velocity field for NGC 4365. Velocities are computed from spectra averaged over Voronoi cells, and are 
corrected for the hs Gauss-Hermite term inside R = I'll. The color bar indicates the velocity scale relative to systemic. Overplotted 
"spider web" shows the polar grid on to which the data are rebinned for model fitting. (6) Rebinned velocity field, antisymmetrized by 
folding and averaging points on opposite sides of the centre. White points mark the bin centres. Colors show a continuous velocity field 
interpolated from the rebinned data. The whited-out central region is excluded from the fit. Dotted lines show V-band isophotes. This 
and all similar figures have North up and East to the left. 



where f is a scaled angular variable. On a spherical shell, t 
is related to the polar angle 9 by 



t = 



< sin^'^VT, 



> sm 



Vr. 



(1) 



With this definition, < t < 1 is the Z tube segment of 
the boundary and 1 < i < 2 is the X tube segment. We 
take v*{t) to be composed of two piecewise-linear sections, 
described by a Z-tube/X-tube contrast factor C and two 
constants, kx and kz, that describe the variation away from 
the symmetry planes (Figure 0d). The values of C, kx, and 
kz used in the model grid are given in the Appendix. 

Projection of the model at each radius is done using a 
quasi-local approximation. Isophotes are assumed to have 
the shape of the projected iso-luminosity-density surfaces at 
the mean deprojected radius; this assumption is accurate in 
the absence of large shape gradients. The velocity field is 
projected assuming that the three-dimensional flow is sim- 
ilar at radii that contribute significantly to the projection 
integral (i.e., close to the tangent point), scaling with radius 
as r~ , and that the luminosity density at these radii scales 
as r~'°. In practice we have found in all previous applica- 
tions that the results are very insensitive to the values of k 
and £, and so in this paper we have simply adopted k = 3 
and £^0. 

Finally, we assume that the iso-mass-density and iso- 
luminosity-density surfaces have the same shape at each ra- 
dius. This is not as strong as assuming that mass follows 
light; our assumption is only that M/L can be expressed as 
a function of local density. But this is an unimportant dis- 
tinction for this application to NGC 4365, since all of the 
data come from inside the effective radius where the dark 
matter contribution is probably negligible. 



2.2 Data 

NGC 4365 was observed with SAURON on the 4.2 m 
William Herschel Telescope in 2000 March. The observations 
and initial data reduction are described by IPavies et alJ 
12001) and iBacon et all i200lf) . More recently, the data 
have been re-redu ced in a homogeneous way for the entire 
SAURON sample jEmsellem et al.j2004|) . This re-reduction 
includes a new spatial binning based on an adaptive Voronoi 
tessellation guaranteeing a mi nimum signal-to-noise ratio of 
60 per bin JCappellari fc Copin 2003) . The data set used 
in this paper is a penultimate version which has a slightly 
lower signal-to-noise ratio per bin, but which is virtually 
indistinguishable from the final version. 

From the re-reduced data set we take the V, a, and 
hs parameters describing the Gauss-Hermite fit to the line- 
of-sight velocity distribution (LOSVD) as derived by the 
Fourier Correlation Quotient (FCQ) method, along with 
their errors. For radii R < 7'.'2, where /13 is significantly 
nonzero, we correct the mean velocity for the skewness of th e 
LOSVD l|van der Marel fc Franx|[l993l : IStatler et alJll996ll . 
These corrections are significant, and reduce the core ro- 
tation by nearly 50%. The corrected mean VF is shown in 
Figure|5Ji, at the full resolution of the SAURON instrument. 

We re-bin the VF on to a coarser polar grid for com- 
parison with the models. This is purely a computational 
expedient, since CPU time and storage requirements scale, 
respectively, with the number of angular and radial bins. 
We are not losing information because there is no signifi- 
cant structure in the VF at scales smaller than the grid. We 
set the centres of the radial bins at 2'.'4, 4'.'0, 6'.'0, 8'.'4, ll'.'2, 
14'.'8, 19'.'2, and 25'.'2. The central I'.'Q is excluded from the 
fit because the VF projection is very inaccurate when the ro- 
tation curve is steeply rising. Each radial bin is divided into 
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Figure 3. Isophota l major-axis position angle (top) and ellipticity (bottom) plotted against mean radius. Data from lForbesI 11994^ and 
iBender et alj il98a) are shown as squares and diamonds, respectively. Solid and dashed lines show means and Itr errors, resectively, over 
the radial bins used for model fitting. Dotted lines indicate fits to the gradients over each bin, which contribute to the errors as explained 
in §1^ 



between 10 and 24 approximately square angular segments, 
which contain between 4 and 12 Voronoi nodes. The grid 
is shown overplotted in Figure |5Ji. We also antisymmetrize 
the VF by folding and averaging points on opposite sides 
of the centre. The symmetric component, i.e., the difference 
between the rebinned data and the folded version, has an 
RMS amplitude of 3.3kms~^, which is consistent with zero 
given the observational errors. We use the symmetrized VF 
to compare with the models, since the latter are by con- 
struction point-symmetric. The antisymmetrized, rebinned 
VF is shown in Figure I^Jj. 

Profiles of isophotal ellipticity and position angle 
are ta k en fr om the HST/PC F555W photometry of 
Forbej ll99^ and ground-based '[/-band photometry of 
Bender et alj il988l) . The profiles are plotted in Figure |21 
Combining the datasets is slightly complicated by the fact 
that neither set of profiles was published with error bars; 
the error bars in the figure represent a guess as to the qual- 
ity of the data. We discard the ground-based ellipticity data 
for R < 9.6", because it differs significantly from the HST 
data in the sense expected from seeing effects. We average 
the remaining data over each radial bin, and also fit for a 
gradient across the bin. The adopted errors represent the 
formal error in the mean added in quadrature with 19% of 
the systematic variation across the bin. The dotted black 
lines in Figure |2J3 show the isophotes, plotted at mean radii 
matching the radial bin centres. 



2.3 Fitting 

Model fitting is done using Bayesian methods, which yield 
probability distributions for the axis ratios at each radius 



and the orientation of the system. The procedure is concep- 
tually straightforward: for each computed model, the likeli- 
hood of obtaining the observed velocity field and isophotal 
ellipticity and position angle profiles is calculated from the 
model prediction and the known measurement errors. This 
is repeated for ~ 10* models covering the space of shape, 
orientation, and internal dynamical parameters, creating a 
multidimensional likelihood function. The likelihood is mul- 
tiplied by a model for the parent distribution of these pa- 
rameters (the "prior"), and integrated over the parameters 
that one is not trying to constrain. The normalized result 
is the so-called "posterior" probability distribution for the 
parameters of interest. 

In practice, the procedure is slightly complicated by 
technicalities needed to handle the 42-dimensional param- 
eter space. The likelihoods for each radial bin are stored on 
a 20 X 20 grid in triaxiality T and axis ratio c/a, and an angu- 
lar grid with approximately square bins of ~ 10° resolution. 
The contribution to the likelihood from the VF and elliptic- 
ity data can be calculated independently at each radius, re- 
ducing storage demands considerably. We find it desirable at 
this stage to introduce a penalty for configurations that re- 
quire unrealistically high internal velocities because the line 
of sight is close to the rotation axis (see the Appendix for a 
full explanation) . Turning the likelihoods into the shape pro- 
file of the galaxy formally involves the construction of a joint 
distribution in 16 variables. To make this more manageable, 
we compute the projections of this distribution into each 
of the shape planes, and use these "marginal distributions" 
as the basic description of the shape profile. The procedure 
for combining the likelihood fits at each radius and including 
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Figure 4. Intrinsic shape profile and orientation of NGC 4365, based on tlie "maximal ignorance" unweighted average of all models 
and a flat parent shape distribution. Rectangular panels show probability densities, in the space of triaxiality T and short-to-long axis 
ratio c/a, for each radial zone; radii are indicated above each panel. Round figure at lower right shows probability density for the line 
of sight over a hemisphere, in Lambert equal-area polar projection; centre, right and left edges, and top and bottom edges of the figure 
correspond to views down the short, middle, and long axes of the galaxy, respectively. The direction of the angular momentum vector 
corresponds to a point along the upper half of a vertical line through the centre. The greyscale is logarithmic in all panels, and white 
contours enclose 68% and 95% of the total probability. 



th e isophotal PA data is explained in detail in the Appendix 
of lStatler et~all JlQQgtl . 

Construction of the prior involves a few subtle points 
which are described in detail in the Appendix of this paper. 
In brief, the prior distribution is assumed to separate into 
four factors: 



F^ 



= -^F^^r[T^, {c/a),]F,or[T^ - T,, (c/a), - {c/a)j 



4-K 



t[Ci,kxx,kzi;Tt, {c/a)i] 



(2) 



The leading factor is the angular part, which we take to be 
isotropic over the viewing sphere. The second is the par- 
ent shape distribution; in principle this is a model for the 
shape distribution of all ellipticals, although for modeling 
a single system it is arguably better to use a flat distribu- 
tion. We will show results for both a flat distribution and 
for t he parent distributio n derived from published long-slit 
data JBak fc Statleil200CI) . The third part describes correla- 
tions in the joint shape distribution among the radial bins. 
These correlations reflect the fact that we expect the intrin- 



sic shape to change continuously with radius. The last factor 
describes the distribution of the internal parameters, C, kx, 
and kz ■ In the "maximal ignorance" estimate we assume no 
correlations of kinematics with radius. In other words, we 
explicitly allow for the possibility that there can be kinemat- 
ically different populations dominating at different radii, bu 
we make no attempt to distinguish these populations, be- 
yond their effects on the local mean streaming motions. In 
our standard model grid we also include a subset of models 
in which C is correlated with shape at each radius, but these 
models give results consistent with the others. 



3 RESULTS 

3.1 Intrinsic Shape Profile 

The basic results are shown in Figure 01 The rectangular 
panels show the posterior probability distributions for the 
shape of the isodensity surfaces in the eight radial bins, 
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Figure 5. Same as Fig.|^ but using the parent shape distribution from lBak fc Statlerl 12001 



based on an unweighted average of all computed models. 
Each panel depicts the (T, c/a) plane and is plotted so that 
oblate spheroids {T = 0) he along the right edge and pro- 
late spheroids (T = 1) lie along the left edge. Spheres occupy 
the entire top margin at c/a — 1, and the bottom margin, 
at c/a — 0.35, corresponds to the observed upper limit on E 
galaxy ellipticities. The white curves indicate the 68% and 
95% highest posterior density (HPD) regions, i.e., the con- 
tours enclosing 68% and 95% of the total probability, and 
may be thought of as la and 2a error regions. 

It is clear that NGC 4365 is an almost maximally tri- 
axial system. The preferred triaxialities are in the vicinity 
of T « 0.4-0.5, and axisymmetric and near-axisymmetric 
shapes are strongly ruled out, except perhaps in the inner- 
most bin. There is an indication of an outward triaxiality 
gradient. The galaxy is also quite a bit flatter than it ap- 
pears; photometry alone requires only that c/a ^ 0.75. 

The results in Figure |1] are for a flat prior shape distri- 
bution, which is appropriate for one system modeled in iso- 
lation. But this part of the prior actually represents the par- 
ent shape distribution for the full population of ellipticals, 
whi ch can be estimated from publis hed long-slit kinemat- 
ics. | Bak fc Statlen i2000|) model the iDavies fc Birkinshawl 
il988ll sample of 13 ellipticals, using methods very close to 
those of this paper, and derive several different parent dis- 



tributions based on different assumptions for the dynami- 
cal prior. Their "maximal ignorance" result (their Fig. 2a) 
corresponds most closely to the dynamical prior used here. 
This distribution has broad peaks centred on the oblate and 
prolate limits, approximately at (T, c/a) = (0, 0.68) and 
(1, 0.82), with a broad ridge between them. Spheres and very 
flat, prolate-triaxial shapes have very low amplitude. 

The shape proflle of NGC 4365 using the Bak & Statlej 
(2000) parent distribution is shown in Figure |S| The results 
are similar to those using the flat parent distribution, dif- 
fering primarily in that the extended tail toward large T 
and small c/a is lacking. It is worth noticing that the single- 
galaxy shape constraint is improved when information from 
other objects is included in the form of the parent distribu- 
tion. One can expect that the greatly superior parent dis- 
tribution that will be derivable from the SAURON sample 
will permit even more precise single-galaxy shapes to be de- 
termined. 

One dimensional triaxiality proflles can be straightfor- 
wardly calculated by integrating the distributions in Figures 
2] and 13 over c/a. The results are shown in Figure |S| Error 
bars show the 68% HPD intervals, which are nearly centred 
on the expectation values. The systematic effect of changing 
the parent shape distribution is smaller than the statistical 
errors. There is an indication of a weak triaxiality gradient. 
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Figure 6. Radial triaxiality profiles derived from Figures El andlSl 
for the fiat and Bak & St atler ( 2000) parent distributions. Points 
with error bars show expectation values and 68% HPD intervals 
(Icr error bars); horizontal bars indicate the widths of the bins. 
Solid lines below denote 99%-confidence lower limits on T. 
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but constant T cannot be excluded. Possibly of greater in- 
terest for constraining the effects of chaos are lower limits on 
T; the continuous curves in the figure show 99%-confidence 
limits, which are significantly non-zero except possibly in 
the innermost bin. 



3.2 Orientation 

The round panels in the lower right of Figures 01 and |S] show 
the posterior distributions for the line of sight, over one 
hemisphere.^ These are the orthogonal projections of the 
full posterior probability, integrated over shape, in contrast 
to the shape distributions which are integrated over orien- 
tation. 

The line of sight is constrained to two fairly narrow 
strips, in the quadrant containing the angular momentum 
vector. That is, the true rotation axis in the outer part of the 
galaxy is most likely pointing in the general direction of the 
observer, as opposed to lying in the plane of the sky. Lines 
of sight ~ 30° from the long axis are preferred, and those 
near the short axis are strongly excluded. These constraints 
are only slightly afi'ected by the choice of the parent shape 
distribution. 

The two strips in the orientation figures correspond to 
views on opposite sides of the {x, z) plane, producing the 
same velocity field. Mirror-image points on the left and right 
sides of the figure differ only photometrically: for a given 
triaxiality gradient the isophotes twist in the opposite di- 
rection. If there were no observed isophotal twist, the figure 
would be left-right symmetric. The small asymmetry in the 
figure reflects the small (~ 2°) observed twist over the fitted 
region. 



■^ Views from antipodal points on the sphere are mirror images 
of each other. The likelihoods recorded in each angular bin are 
averages of the two antipodal views. The ambiguity is trivially 
resolved a posteriori. 



Figure 7. Lambert equal-area projections of half of the view- 
ing hemisphere, showing (top) T and (bottom) c/a values of the 
best-fitting optimized model at each orientation (black contours) . 
Overplotted in red are contours of Axponi ^^ levels of 2.3, 6.17, 
and 11.8. Models outside the outermost red contour are formally 
excluded at > 99% confidence. Dotted lines show an angular grid 
at 10° spacing for reference. Letters in the lower panel indicate 
the 3 models shown in detail in Fig. |S] 



3.3 Goodness of Fit; Ax^ Confidence Limits 

A disadvantage of the Bayesian approach is that a poste- 
rior probability can be calculated even if the fit to the data 
is poor. It is therefore necessary to check separately that 
the best models are actually adequate fits in a x^ sense. 
We couple our modeling code to a standard multidimen- 
sional optimizer and calculate, for each angular bin centre, 
the combination of shape and dynamical parameters that 
minimizes x^ at that orientation. This allows us also to ver- 
ify whether the HPD regions agree with confidence limits 
determined from y^ contours. 

We have 81 data points (66 velocities, 8 ellipticities, and 
7 PA differences), and, at fixed orientation, 40 parameters 
(T, c/a, C, kx, and k^ in each radial bin). We minimize 
a x^ statistic that includes the same penalties (for shape 
differences between adjacent bins and large internal v/a) 
that are used in the Bayesian treatment. This adds 15 con- 
straints, making i/pen = 56 degrees of freedom. The best 
optimized model has Xpcn/'^ = 1-14. This is an acceptable 
fit; the probability of a larger Xpon occurring by chance is 
0.22. We also calculate, but do not optimize, the unpenalized 
xSn nieasuring only the fit to the data. The model optimized 
on Xpon has Xun/^'un = 1.41, indicating an adequate but not 
spectacular fit. 

In Figure|7]we plot in black the mean shape (unweighted 
average over the 8 radial bins) of the best optimized model at 
each orientation, for the half of the viewing hemisphere con- 
taining the viable models. These contours indicate only the 
best model for each line of sight, not whether these models 
are statistically good fits. Overplotted in red are contours 
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Figure 8. Data (top left) and models in the observational domain. Top row: isophotes and velocity fields, plotted as in Fig. El White 
points mark bin centres, where the fit is performed. Bottom row: ray-traced images of the mean shapes, as oriented in the sky. The 
shapes are rendered as polished metallic ellipsoids, symmetrically illuminated by lights to the right and left behind the observer. Red, 
green, and blue rods indicate the long, middle, and short axes. Models a, b, c correspond to those indicated in Fig.|7| with (a) good, (b) 
marginal, and (c) poor fits. 



of Axpcn, at levels of 2.3, 6.17, and 11.8 above the mini- 
mum. Formally these correspond to 68%, 95%, and 99.7% 
confidence intervals for two degrees of freedom, so lines of 
sight outside the last red contour are strongly ruled out. 
The models inside the red contours, which are not ruled 
out, have triaxialities T between 0.32 and 0.72, and fiatten- 
ings c/a between 0.40 and 0.68. This is consistent with the 
results in Figure |1] which include many more models than 
just the best ones on each line of sight. Furthermore, the 
correspondence between the red Axpen contours in FigureQ 
and the HPD contours in the bottom right of Figure 2] is 
qualitatively good; they clearly identify the same preferred 
hues of sight. The agreement is not exact, however. We at- 
tribute this to the difference between Ax^ values computed, 
on the one hand, along an optimised x^ valley floor in the 
space of the non-orientation parameters and, on the other, 
averaged over the whole valley at quite coarse resolution. 

Notice that, on the viewing hemisphere, the x^ valley 
is nearly parallel to the contours of optimized c/a, but cuts 
across the contours of optimized T. The fact that triaxialities 
~ 0.45 are optimal over the broadest range of orientations is 
consistent with those values of T being most probable (Figs. 
13] |S] andlSJ. One should not place too much emphasis on the 
exact position of the x^ minimum. Visual inspection of the 
acceptable models shows that different models fit different 
parts of the data better than others, and an alternative re- 
binning of the VF might easily have moved the x^ minimum 
by ~ 10° along the valley. 

The blue letters in Figure |7| identify three representa- 
tive models which we show in detail in Figure |5] Models a, 
b, and c correspond to good, marginal, and bad fits respec- 
tively, having Xun/^^ values of 1.4, 1.7, and 2.8, and Xpon/'^ 



values of 1.1, 1.3, and 2.3. The top row of Figure |5] shows 
the projected isophotes and velocity field of each model, 
with the observed data from Figure |5j3 reproduced at top 
left. The bottom row shows a rendered-surface representa- 
tion of an ellipsoid with the mean shape of each model, at 
the correct orientation. Clearly the acceptable models are 
able to reproduce the detailed kinematic and photometric 
structure of the galaxy. One can see from the figure that the 
biggest difficulty is simultaneously fitting the rotation am- 
plitude and the steep velocity gradient near PAs 100° and 
280°. The unacceptable model c does a particularly poor job 
of reproducing this gradient, and its kinematically distinct 
core region is too small. 



3.4 Orbital Structure 

Because we have allowed the orbit populations, by way of 
the velocity boundary condition, to vary freely with radius, 
the inferred shape profile and orientation of the system are 
essentially free of any assumptions about dynamical sub- 
systems in the galaxy or its formation history. But, by the 
same token, we are unable to constrain the details of the 
stellar populations with these coarse models. That would be 
a job for much more sophisticated modeling, making use, 
at least, of the full line-of-sight velocity distributions, and 
possibly line-strength indices as well. None the less, we can 
draw some conclusions concerning what orbits dominate the 
net streaming motions in the inner and outer parts of the 
galaxy. 

To explore the dynamics of the kinematically distinct 
core, we compare three subsets of models in which the in- 
nermost zone is restricted to certain values of the contrast 
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Figure 9. Intrinsic axis ratios of merger remnants and dark haloes from the literature. Each panel shows mean shapes or local shapes 
near the half-mass radius, as reported by the authors. Left: spiral-spiral merger remnants, without (Hernquist 1992, 1993, Barnes 
1992) and with (Barnes & Hernquist 1996) gas dissipation. Center: Remnants of dissipationless group mergers (Weil & Hernquist 1996 
as reanalysed by Statler et al. 2001, Zhang et al. 2002). Right: Distribution of zero-redshift halo shapes in a dissipationless ACDM 
cosmological simulation (Jing & Suto 2002) {contours; spaced by factor of 2 in density) and individual shapes of luminous objects in 
dissipational simulations (Sugarman et al. 2000, Meza et al. 2003). 



C: a Z-tube dominated configuration (C = oo, 2), an X-tube 
dominated configuration (C = 0, 0.0625), and a neutral con- 
figuration (C = 0.25, 0.5)'* We compute the maximum of 
the joint likelihood L{T, c/a, cf>, 6) in each radial bin, inte- 
grated over dynamical parameters but not over orientation. 
We find the Z-dominated configuration to have the highest 
peak likelihoods, with those in the neutral and X-dominated 
configurations smaller by factors of 5 and 3 x 10*, respec- 
tively. Insisting that the core rotation is entirely due to Z- 
tubes barely alters the shape profile. Axisymmetry in the 
core is still excluded at a > 2cr level. Photo metrically, the 
isophotes between 1" a nd 3" are slightly discy (ICarollo et alJ 
ll997l : lRest et al]|200J) . and the large /13 values in the same 
region are suggestive of discy kinematics. But a stellar disc 
embedded in a larger triaxial system will not, in general, be 
axisymmetric. 

We can also investigate the dynamics of the region out- 
side the core, by following the same procedure and restrict- 
ing C in the outer 4 zones. Here we find the situation re- 
versed: the X-dominated configuration has the largest peak 
likelihoods, with the neutral and Z-dominated configura- 
tions down by factors of 4 and lO'' , respectively. 

We conclude, therefore, that the mean motions are most 
likely dominated by long-axis tubes in the main body of 
the galaxy and by short-axis tubes in the kinematically dis- 
tinct core; this is essentially the configuration advocated by 
[Surma fc Bender (1995) for NGC 4365, and by Franx ot ai.i 
il99ll) for the similar system NGC 4406. It is virtually cer- 
tain that the main body is not short-axis tube dominated, 
and that the core is not long-axis tube dominated. But con- 
figurations in which both tube families contribute signifi- 
cantly all through the galaxy cannot be ruled out. In fact, 
all of the models in Figure |5] require Z tubes in the main 
body to move the apparent rotation axis away from the po- 
sition angle of the projected long axis. 



^ Z tubes still dominate the total velocity field when C 
can be seen in Fig. 



1, as 



4 DISCUSSION 

4.1 The Importance of Quantitative 
Measurements of Galaxy Shape 

The distribution of three-dimensional shapes is an under- 
utilized, but potentially powerful diagnostic of galaxy for- 
mation physics. Numerical simulations of structure forma- 
tion and galaxy mergers can now predict the axis ratios of 
the zero-redshift haloes or final merger remnants, and the 
results depend on the processes involved. Figure|n|shows rep- 
resentative results from the literature. Pair mergers of pure 
discs produce very flat, very triaxial remnants (Hornauisa 
[1992). Including dense bulges or moderate gas dissipation 
tends to mak e the final systems rounder and less prolat e 
iBarnejll993 : lHernauisjll993 : iBarnes fc HernauistlllQQ^) . 
in the former case because the merging systems are them- 
selves rounder, and in the latter because of the scattering 
of low angular momentum orbits by the central mass ac- 
cumulation. Significant dissipation in gas-rich mergers can 
fiatten the final system if the gas settles into a disk and 
subsequently forms stars (Hernauist & Barnes 1991; Barnea 
2002; Lamb fc Hearn 2003) . Pure dark matter haloes formed 
in cosmological sim ulations are also very flat and prolate 
JJing fc Satdl2002) : as in disc mergers, including some dissi- 
pation tends to m ake the final systems rounder, though sti ll 
strongly triaxial ISuga rman et al. ' 2000; M eza et al.l 12003) . 
Most prominently, mergers of groups should produce highly 
prolate or oblate remnants, dependi ng on whether the initial 
system is co llapsing or virialized fWe il fc Hernauis1ill996l : 
iGariio et aTl llQQT: Zhang ct al. 2002). Clearly there is dis- 
criminating power in the shape distribution. 

The SAURON instrument has produced stel- 
lar kinematic maps for two dozen elliptical galaxies 
JEmsellem et al.] |200J). This paper constitutes the first 
step in using these data to determine the shape distribution. 
The results of the preceding section illustrate the power of 
two-dimensional data to constrain the shapes of individual 
systems. The axis ratios of NGC 4365 are determined 
to within approximately 0.1, at the Icr level. We have 
also modeled the SAURON VF for NGC 3379, and find 
that we can constrain the shape to comparable precision. 
The NGC 3379 results are v ery simi l ar to those obtained 
from 7 long-slit profiles by IStatlerl i200 J) . This is not 
surprising, since the VF is very simple and symmetric. 
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and the SAURON map does not contain much additional 
information beyond what could be interpolated from the 
long-slit data. One may justifiably anticipate that it will be 
possible to determine the shapes and orientations of other 
objects in the SAURON sample to comparable precision 
using these methods, to estimate the paren t shape distri- 
bution for the sample with the approach of iBak fc Statlerl 
'200(f), and to study the dependence of this distribution on 
galaxy luminosity or environment. These topics will be the 
focus of future papers. 



4.2 Effects of Chaos; Limits on Central Black 
Hole Mass 

NGC 4 365 is an old system. From the metal and Balmer line 
indices iDavies et al.l ll2001f) infer an age ^ 12 Gyr for the 
bulk of the stellar population, with at most a few percent 
admixture of a younger (5 Gyr) population in the centre. Ab- 
solute s tellar ages are unc ertain, however, and recent photo- 
metr ic JPuzia et al.l 12003) and spectroscopic JLarsen et alJ 
l2003ll studies of globular clusters indicate the presence of an 
intermediate-age population. A minor merger 2-5 Gyr ago 
is not unlikely, although a significant event would probably 
have left photometric t races at large radii, which are not 
seen JDavies et alJl200]J) . 

None the less, any age in the range of several Gyr im- 
plies that the system is dynamically very old. We can define 
a local relative dynamical age rd{r) as the ratio of the as- 
sembly age t* to the crossing time at projected radius r. 
We approximate the latter as r/[v^{r) + (T^(r)]^", where 
Vm{r) is the maximum mean velocity at r. The result can 
be represented extremely well by the relation 
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thus the entire SAURON field is hundreds of dynamical 
times old even if the galaxy was assembled as recently as 
2 Gyr ago. 

The long-term survival of triaxiality in systems with 
central density cusps or black holes h as been questioned 
many times JLake fc NormarJ Il983t iGerhard fc BinnevI 
Il985l : iMerritd Il999l . and references therein). Central mass 
concentrations can render a significant fraction of the box 
orbit phase space chaotic. Since box orbits are essential for 
maintaining triaxiality, it is believed that growth of a cusp 
or black hole can force a triaxial system to evolve toward 
axisymmetry or sphericity. That this has clearly not hap- 
pened in NGC 4365 allows us to limit the maximum mass 
of a central black hole. 

Calculating a precise upper limit is not yet possible. 
Surprisingly few numerical studies of chaos-driven mor- 
phological evolution have been performed, and these are 
jiot entirely appropriate to the properties of NGC 4365. 
iMerritt fc Ouinlcml lll99a) start with a prolate-triaxial {T « 
0.6, c/a « 0.5) equilibrium system with a fiat, constant- 
density core, and grow central dark masses amounting to 
0.3%, 1%, and 3% of the galaxy mass. They see a large, rapid 
decrease in T that effectively axisymmetrizes the entire sys- 
tem for mass ratios > 1%. Below this mass the change is less 
drastic, but still pronounced. They also find that the shape 
evolution continues even a fter the central object reaches its 
final mass. lSellwoodI l|2003) argues that this late evolution is 



a numerical artefact, but confirms Merritt fc Quinlan's ba- 
sic result that mass ratios larger than 1% globally destroy 
triaxiality. iHoUev-Bockelmann et aU (|2002) perform a simi- 
lar experiment in which a 1%-mass black hole is grown in a 
somewhat rounder (T = 0.54, c = 0.70) triaxial system with 
an initial r~^ density cusp. Unlike the flat-core systems, the 
cuspy system does not experience a global destruction of 
triaxiality. Its inner regions tend toward a spherical shape 
at approximately constant T, while the outer half of the 
system remains flattened and triaxial. The black hole cre- 
ates a spherical region around itself containing several times 
its own mass in stars. The interpretation of these simula- 
tions is further complicated by the existence of flattened 
triaxial equilibria cont aining both cusps and central masses. 
iPoon fc Merritd (120021) construct equilibrium models for tri- 



axial (T = c/a = 0.5) nuclei with r~ and r~ cusps by 
Schwarzschild's method. They confirm using A'^-body simu- 
lations that the shapes, even down to radii enclosing only a 
few times the black hole mass, are stable for ~ 10 crossing 
times. It is not known whether they would persist over tens 
of thousands of crossing times, as NGC 4365 would require. 
For whatever reason, these equilibria seem not to be found 
by the simulations that begin with no black hole. Evidently 
they do not have strong basins of attraction around them 
when the central mass grows dynamically. 

W e accept the basic result fr om [ Merritt fc QuinlM 
19981) . iHollev-Bockelmann et ai] ^2002^ ■ and Sellwoo" 



20021) ■ that a black hole > 1% of the system mass will 



drive a rapid evolution of a region whose size is a significant 
fraction of the effective radius toward either axisymmetry 
or sphericity. This can be decisively ruled out in NGC 4365. 
To estimate the system mass, we use t he total 1^-band mag- 
nitud e, V^ = 9.5, given in the RC 3 ide Vaucou leurs et all 
Il99ll) . a distance modulus of 31.55 JTonrv et alJl2001^ . and 
a mass-to-light ratio MjLv = 6 obtained by fitting an 
isotropic oblate Jeans model to the full SA URON data set. 
This M/ Lv is consistent with the results of lCebhardt et alJ 
f2OO30, who find a mean MjLy = 5.3 for 11 ellipticals, 
with a dispersion of 2.1. We obtain M = 3.3 x 10^^ M© 
for NGC 4365. This implies that a black hole of mass 
Mbh > 3 X 10^ Mq would either globally axisymmetrize 
the galaxy or at least render the inner 10^^° Mq spherical. 
By integrating the 7 = 0.13 Nuker mo del fit to the sur- 
face brightness profile IIRest et al. l 2001), we find that this 
mass is enclosed in projection by the isophote with a mean 
radius of 3". Although we see indications of declining triax- 
iality at these radii, there is no significant rounding of the 
isophotes, which would be seen if the inner regions were be- 
coming spherical. We conclude that the long-lived triaxiality 
in NGC 4365 rules out black holes of Mbh > 3 x 10^ Mq. 

There is, as yet, no direct evidence for a black hole in 
NGC 4365. But for a moan dispersion of 260kins ~^, the 
MgH-o" relation (Gobhardt ot al. 2000; Fcrrarcsc fc Merritd 
I2OOOI: iTremaine et alJ l2002l) would predict Mbh ~ 4 x 
10* Mq. This corresponds to ~ 0.1% of the system mass, 
and probably would not drive a global shape change. If 
the single simulation of lHoUev-Bockelmann et al.l i200 2') can 
be extrapolated to lower masses, then a black hole of this 
mass might be expected to affect the isophotes in the inner 
0'.'9. In the simulation, an r~^ density cusp grows within 
the classical radius of influence of the black hole, which, 
for NGC 4365's central dispersion of 270kms~^, would 
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subtend (/.'l. HST /WFP C2 surface photometry (iRest et alJ 
I2001I : ICarollo et al.lll9iMI does show a gradual rounding of 
the isophotes interior to l", reaching a minimum ellipticity 
~ 0.08 at about O'.'IG. The central photometric structure is 
complicated, and may be affected by dust. There is no steep- 
ening of the 7?- band brightness profile within 0.1", but there 
does seem to be a central, unresolved source t hat is slightly 
(^ .02 mag) bluer than its surroundings llCarollo et alJ 

Thus, while the persistent triajciality of NGC 4365 ex- 
cludes large black holes, a more modest object, on the scale 
predicted by the Mbh-o" relation, is probably small enough 
not to globally alter the shape of the galaxy and can be 
consistent with the data. 



4.3 Figure Rotation 

We have shown in i; l3.4l that net streaming motions in long- 
axis tubes contributes significantly to, and probably domi- 
nates, the observed velocity field. The direct and retrograde 
branches of these orbits generate the same spatial density 
distribution if there is no rotation of the potential. How- 
ever, if the galaxy is tumbling about the short axis, these 
orbits acquire a tilt with respect to the symmetry plane. 
This tilt increases for larger orbits, and is in the opposite 
direction for two senses of circulation. Since there must be 
an imbalance in the populations of the direct and retrograde 
branches, fi gure rotation would induce an intrinsic twist in 
the galaxy llSchwarzschildll982r . The small observed photo- 
metric twist allows us to limit the tumbling rate. 

The tilt of the periodic planar orbits out of the {y, z) 
plane c aused by figure rot a tion a bout the z axis wa s exam - 
ined bv Ivan Albada et al.l (Il982h and iHeisler et al.l Jl982ll . 
These so-called "anomalous" orbits are the parents of the 
long-axis tubes ^gchwarzschild 1982), which we will as- 
sume acquire roughly the same tilt as the parent of the 
same energy. This assumption is by no means firmly estab- 
lished; the behavior of general orbits in tumbling figures has 
been studied only superficially (Salow 1998, unpublished). 
Ivan Albada et alJ ( j.1982) show that the tilt angle a can be 
approximated by 

2 ^ 

1 — {k/ujY + {il/ijj) oj 

where uj and k are the frequencies of the periodic motion in 
the plane and of oscillations out of the plane, respectively, 
in the absence of figure rotation, and Q, is the tumbling fre- 
quency. For slow tumbling, SI can be ignored in the first 
factor. The ratio k/uj is a function of the axis ratios, and 
should be smaller than 1 for orbits around the long axis. 
For a triaj cial modified Hubble mo del with T = 0.8 and 
c/a = 0.5, Ivan Albada et al.l il982fl calculate k/uj ~ 0.65. 
Our models for NGC 4365 suggest shapes not quite as ex- 
treme as this. Let us suppose that k/u) ~ 0.8, which is ap- 
propriate for a prolate logarithmic potential with a density 
axis ratio of roughly 0.6. For NGC 4365 we approximate uu 
by 2^'^a/r, with a ~ 220kms~^ the one- dimensional dis- 
persion at i? ~ 20" , and find that 
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where Pn = 2n/Q is the tumbling period. 



It is interesting that the tilt is predicted to be linear in 
r. NGC 4365, in fact, has a linear isophotal twist of 0.09°/" 
for 2" < R < 60". How the anomalous orbit tilt translates 
into an observed isophotal twist is extremely complicated, 
depending on the exact orbit populations, axis ratios, and 
viewing geometry, and is beyond our ability to predict in 
detail. However, if we imagine that the observed twist and 
the anomalous orbit tilt are of the same order of magnitude, 
then the observed twist could be caused by figure rotation 
if the tumbling period Pq ~ 5 Gyr. Alternatively, we can 
say that Pq ^ 500 Myr is probably ruled out by the absence 
of a larger twist. This tumbling rate would place corotation 
at about 8 effective radii assuming a standard isothermal 
dark halo. It seems, therefore, safe to conclude that figure 
rotation is dynamically unimportant in the observable part 
of the galaxy. 

If the isophotal twist is indeed due to figure rotation, 
then it is worth asking in wha t direction the figure is ro- 
tating. Ivan Albada et al.l (119821) performed this exercise for 
Centaurus A, in an effort to explain the warp of the dust 
lane. They predicted that the southwest side of the figure 
should be approaching and the northwest receding; unfor- 
tunately this was later found to be opposite to the sense 
of rotation of the stars. For NGC 4365, the required sense 
of figure rotation is northeast approaching and southwest 
receding for models viewed as in Figure |H^ andjSJD. Short- 
axis tubes are needed in these models, to displace the zero- 
velocity contour counter-clockwise on the sky. The sense of 
rotation of the short-axis tubes is therefore in the same di- 
rection as the figure tumbling, but opposite to the orbits of 
the same family that dominate the decoupled core. 

4.4 Implications for Other Modeling Approaches 

Although the VF fitting method can constrain the axis ra- 
tios and the orientation of the mass distribution, it does not 
make use of all of the available kinematic data, and conse- 
quently cannot truly constrain the full density profile. To do 
this, more sophisticated approaches are needed. 

Jeans models have a long and distinguished history in 
galax y dynamics. These models are built on the Jeans equa- 
tions iJeansll919l : IChandrasekhaJl960l : lBinnev fc Tremainj 
Il987l) . which relate the second moments of the velocity 
distribution to the density distribution and the potential. 
Second-moment models are intuitively appealing, sin ce lumi- 
nous ellipticals are supported by velocity anisotropy JBinnevI 
Il978l) . But fully anisotropic models are not easily computed 
because, in general, there are more independent second mo- 
ments than there are Jeans equations to constrain them. 
Very recently, an analytic solution to the triaxial anisotropic 
Jeans equations h as been demonstrated f or systems with 
Stackel potentials ivan den Ven et alll2003f) . It may be pos- 
sible to use this result to create simple, approximate second- 
moment models, and to merge the VF fitting and Jeans tech- 
niques (van de Ven et al., in preparation). 

A full understanding of the dynamical structure of 
a galaxy requires a model for the complete phase space 
distribution function. At present the only workable tech- 
nique to build such a model is Schwarzschild's method 
iSchwarzschildl Il979t) . In this technique a finite library of 
test-particle orbits is populated to reproduce a desired mass 
distribution and fit the photometric and kinematic data. The 
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ne west S chwarzschild c ode for triaxial system s is presented 
bv lVerolm c (2003) and IVerolme et al.l J2004) : but unfortu- 
nately, no direct comparison of Schwarzschild and VF fitting 
results for NGC 4365 is yet possible. Schwarzschild models, 
by allowing any population of stars on any orbit, have far 
more freedom than the continuity-equation models used in 
VF fitting. We suspect that a Schwarzschild code may be 
able to find models that fit the data outside our statisti- 
cally allowed region. But we also suspect that these models 
may be very unsmooth in phase space. How a physical re- 
quirement of smoothness is best impo sed on S chwarzschild 
models is a matter of vigorous debate iValluri et al. 2004). 
The SAURON kinematic maps constitute the best data 
to date for revealing the dynamical structure of elliptical 
galaxies. To fully realize the potential of these data will re- 
quire the application of a variety of modeling approaches, 
including Schwarzschild's method, Jeans models, and VF 
fitting. But it is not feasible to apply the most computa- 
tionally expensive method blindly; each approach should be 
informed by the results of coarser methods that can survey 
more territory. This paper provides an initial rough map 
of the parameter space for NGC 4365, as a guide for more 
detailed modeling efforts to follow. 



served photometric twist we infer that the tumbling period 
is ^ 500 Myr, putting corotation at at least 8 effective radii, 
and supporting the common assumption that figure rotation 
is dynamically unimportant in luminous ellipticals. 

NGC 4365 contributes to an emerging picture of the 
role of black holes in giant elliptical galaxies: though om- 
nipresent, black holes are not massive enough to alter the 
global structure of their hosts, which can remain anisotropic 
and triaxial for many hundreds of dynamical times. It is 
still undetermined whether this same picture may apply to 
lower-luminosity systems, however, and it is possible that 
there is a threshold below which central black holes drive 
a secular evolution toward axisymmetry. The SAURON 
project has produced a rich database which will make it 
finally possible to address these and related issues. Deter- 
mining the mass distributions and the dynamical structure 
of early-type galaxies will require extensive modeling using 
Schwarzschild's method, as well as faster but more idealized 
moment-based methods such as we have used. We hope that 
through an application of such complementary approaches, 
it will be possible to develop a more complete understanding 
of the dynamical histories of hot stellar systems. 



5 SUMMARY 

We have modeled the two-dimensional mean velocity field 
and isophotal profiles of the old and kinematically inter- 
esting elliptical galaxy NGC 4365 using the VF fitting ap- 
proach, and constrained the system's intrinsic shape be- 
tween 0.03 and 0.5 effective radii and its orientation in space. 
We find that the galaxy is highly triaxial ((T) ss 0.45) and 
somewhat flatter than it appears {{c/a} « 0.6). Axisymme- 
try or near axisymmetry (T < 0.1) is ruled out at better than 
95% confidence everywhere, and at better than 99% confi- 
dence everywhere except in the inner 3". There is a weak 
indication of an outward triaxiality gradient. The orienta- 
tion of the galaxy is constrained to two fairly narrow strips 
on the viewing hemisphere. In the most probable models 
the system is oriented with its long axis pointing in the gen- 
eral direction of the observer, extending in projection to the 
southwest. We have verified that the best models are statis- 
tically acceptable in a x^ sense, and that Ax^ regions and 
Bayesian HPD regions define essentially the same region of 
permitted parameter space. 

That strong triaxiality has evidently persisted in the 
galaxy for some 12 Gyr implies that any central black hole 
must be smaller than 3 x lO'' -Af0 . A larger central mass 
would have induced chaos sufficient either to globally ax- 
isymmetrize the galaxy or at least make the inner several 
arcseconds spherical. On the other hand, a black hole of 
mass 4 x 10* Mq , as predicted by the A/bh-o" relation, would 
probably not preclude long-lived triaxiality and is consistent 
with the observations. 

We find that there must be a significant contribution 
from long-axis tubes outside the kinematically decoupled 
core, and that the direct and retrograde branches of these 
orbits must be unequally populated so as to produce the 
observed rotation signature. This permits us to constrain 
the rate of figure rotation about the short axis, since fig- 
ure rotation would induce an intrinsic twist. From the ob- 
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APPENDIX A: PRIOR DISTRIBUTIONS AND 
PENALTY FUNCTIONS 

The separation of the prior distribution into factors for the 
orientation, the parent mean shape distribution, shape cor- 
relations, and the dynamical parameter distribution is ex- 
plained in § 12.31 This appendix describes the details of the 
last two factors. 



Al Shape Prior 

The prior shape distribution for 8 radial bins is a joint dis- 
tribution in 16 variables. It may seem least biased to let 
this function be constant everywhere and allow any shape 
at each radius, but this is not the case. The shape param- 
eter space at each radius is two-dimensional, and so a flat 
distribution implies that large differences in shape between 
adjacent radii are likely simply because there is more area in 
a plane far from a given point than close to it. If the obser- 
vations show that the galaxy appears similar at neighboring 
radii — in particular, t hat the isoph otal twist is small — then 
the Bayesian "robot" (|Javnesl2003r is forced to reason as fol- 
lows: neighboring radii probably have very different shapes. 
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yet in projection they appear very similar; so we must be 
looking at the galaxy along one of the special lines of sight 
where big changes of shape are not noticeable. The results 
will then be biased toward lines of sight in the symmetry 
planes. 

From observations of many ellipticals it is already 
known that radial photometric gradients tend to be small, 
and therefore that shape profiles are probably smooth and 
shapes at neighboring radii are correl ated. Expecte d corre- 
lations must be included in the prior (|,Tavnesll2003f) . We do 
this by imposing a Gaussian penalty factor of the form 



^ COT 



n 



exp 



(r, -r,-i)^ + [(e/a), -(e/g),- 
2a| 



,(A1) 



which penalizes differences in T and c/a between neigh- 
boring radii (SSt atlej liooH § 3.2.2). We set the constant 
ffg = 0.126, so as to render a complete change of shape, 
e.g., from T = to T = 1 or c/a = to c/a = 1 over the ob- 
served range of radii, a priori unlikely at the 3a level. Other 
values of ag in the range 0.1 to 0.2 give very similar results. 



A2 Dynamical Prior 

The parent distribution of the dynamical parameters C, 
kx, and kz is defined by the grid of values over which we 
compute models. For the Z-tube/X-tube contrast we take 
C = (0, 0.0625, 0.125, 0.25, 0.5, 1.0, 2.0, oo). The limiting val- 
ues C = and C = oo correspond to zero net streaming in 
Z-tubes and X-tubes, respectively. We also add four more 
cases, where C is assumed to be a function of T and c/a 
according to each of the four "prescription" functions intro- 
duced in equations (11)-(14) of^atler (1994b|). We find that 
the results from the regular grid alone are nearly indistin- 
guishable from those including the prescriptions, indicating 
that explicit correlations between shape and dynamics are 
unimportant. 

The other constants kx and kz determine the "disc- 
like" or "spheroid-like" character of the mean stream- 
ing in each tube family. A lower value indicates more 
disc-like kinematics, meaning that the mean velocity de- 
creases more rapidly away from the symmetry plane. 
We take kx = (0,0.25,0.5,0.75,1.0,1.25,1.5) and kz = 
(0,0.25,0.5,0.75,1.0), based on the behavior of numerical 
and analytic self-consistent triaxial models [see Figure 8 of 
[Statlcr. a994a) 1. 



A3 v/a Penalty 

The magnitude of any velocity field can be scaled up or 
down by reversing the direction of some fraction of the stars 
on tube orbits. Accordingly, the model velocities are deter- 
mined only to within a constant factor and then scaled to fit 
the data. In previous work we have allowed the scaling con- 
stant to vary freely. However, we find that this practice leads 
to biased results when the apparent rotation axis is far from 
a symmetry axis. The reason is that the fitting procedure 
can most easily produce a model with the apparent rotation 
axis along an arbitrary direction merely by placing the line 
of sight very close to the true rotation axis, and then slightly 
inclining the model in the right direction. But in this case 



the low sin i factor may imply deprojected internal velocities 
that are unrealistically large. 

To avoid these unphysical models, we calculate for each 
projected model at each radius the maximum deprojected 
velocity on the shell, Vs, required by the optimal scaling, 
and then penalize the model by a factor exp{—Vs/2a^). We 
take the parameter a to be the observed velocity dispersion 
at that radius, so in other words we are penalizing models 
with internal v/a values greater than 1. This turns out to 
be a fairly weak penalty in nearly all regions of parameter 
space, affecting lines of sight only within about 15° of the 
true rotation axis. 
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